# load dependencies
library(rgdal)
Loading required package: sp
rgdal: version: 1.4-8, (SVN revision 845)
 Geospatial Data Abstraction Library extensions to R successfully loaded
 Loaded GDAL runtime: GDAL 2.4.2, released 2019/06/28
 Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.6/Resources/library/rgdal/gdal
 GDAL binary built with GEOS: FALSE 
 Loaded PROJ.4 runtime: Rel. 5.2.0, September 15th, 2018, [PJ_VERSION: 520]
 Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.6/Resources/library/rgdal/proj
 Linking to sp version: 1.3-2 
library(geojsonio)
package ‘geojsonio’ was built under R version 3.6.2
Attaching package: ‘geojsonio’

The following object is masked from ‘package:base’:

    pretty
library(ggplot2)
Registered S3 method overwritten by 'dplyr':
  method           from
  print.rowwise_df     
library(leaflet)
library(plotly)
package ‘plotly’ was built under R version 3.6.2Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     

Attaching package: ‘plotly’

The following object is masked from ‘package:ggplot2’:

    last_plot

The following object is masked from ‘package:stats’:

    filter

The following object is masked from ‘package:graphics’:

    layout
library(RColorBrewer)
library(dplyr)

Attaching package: ‘dplyr’

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union
library(broom)
library(mapproj)
Loading required package: maps
library(reshape2)
package ‘reshape2’ was built under R version 3.6.2
library(digest)
library(lubridate)

Attaching package: ‘lubridate’

The following object is masked from ‘package:base’:

    date
library(spdplyr)
library(rmapshaper)
package ‘rmapshaper’ was built under R version 3.6.2Registered S3 method overwritten by 'geojsonlint':
  method         from 
  print.location dplyr
library(raster)
package ‘raster’ was built under R version 3.6.2
Attaching package: ‘raster’

The following object is masked from ‘package:dplyr’:

    select

The following object is masked from ‘package:plotly’:

    select
library(mapview)
package ‘mapview’ was built under R version 3.6.2
library(sf)
package ‘sf’ was built under R version 3.6.2Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0
# try this with: https://www.r-graph-gallery.com/168-load-a-shape-file-into-r.html

zipcode_boundaries <- readOGR("../ZIP_CODE_040114/ZIP_CODE_040114.shp",
  layer = "ZIP_CODE_040114", verbose = FALSE)

# try different table of zipcodes:
#alt_zipcode_boundaries <- readOGR("../alt_nyczip/nyc_zip_code_tabulation_areas_polygons.shp",
#  layer = "nyc_zip_code_tabulation_areas_polygons", verbose = FALSE)

# fortify the data to obtain dataframe for use by ggplot2
zipcode_boundaries_fortified <- tidy(zipcode_boundaries, region = "ZIPCODE")

# plot it
ggplot() +
  geom_polygon(data = zipcode_boundaries_fortified,
               aes(x = long, y = lat, group = group),
               fill = "#69b3a2", color="white") +
  theme_void()
# try using info from: https://nceas.github.io/oss-lessons/spatial-data-gis-law/3-mon-intro-gis-in-r.html

class(zipcode_boundaries)

extent(zipcode_boundaries)

plot(zipcode_boundaries, main = "NYC zipcode boundaries")
# convert spatial data frame to geojson
# nyc_json <- geojson_json(zipcode_boundaries)
# save to local file
# geojson_write(nyc_json, file = "./nyc.geojson")

alt_nyc_json <- geojson_json(alt_zipcode_boundaries)
geojson_write(alt_nyc_json, file = "./alt_nyc.geojson")
# create leaflet object
pal <- colorNumeric("viridis", NULL)

leaflet(alt_zipcode_boundaries) %>%
#  addProviderTiles(providers$OpenStreetMap) %>% 
  addTiles() %>% #add basemap
#  addPolygons(stroke = FALSE, smoothFactor = 0.3, fillOpacity = 1,
#    fillColor = ~pal(borough),
#    label = borough) %>%
  addPolylines(color = "#444444", 
               weight = 1, 
               smoothFactor = 0.5, 
               opacity = 1.0)
### run it all again but this time with SF (faster?)

zipcode_sf <- st_read("alt_nyczip/nyc_zip_code_tabulation_areas_polygons.shp")
Reading layer `nyc_zip_code_tabulation_areas_polygons' from data source `/Users/changhoyoon/Documents/Chang_Ho_Work/Master/Coursework/BMI706 data visualisation/viz_proj/COVID-19_NYC_706/alt_nyczip/nyc_zip_code_tabulation_areas_polygons.shp' using driver `ESRI Shapefile'
Simple feature collection with 262 features and 12 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -74.25576 ymin: 40.49584 xmax: -73.6996 ymax: 40.91517
CRS:            4326
altzip_sf <- st_read("../ZIP_CODE_040114/ZIP_CODE_040114.shp")
Reading layer `ZIP_CODE_040114' from data source `/Users/changhoyoon/Documents/Chang_Ho_Work/Master/Coursework/BMI706 data visualisation/viz_proj/ZIP_CODE_040114/ZIP_CODE_040114.shp' using driver `ESRI Shapefile'
Simple feature collection with 263 features and 12 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 913129 ymin: 120020.9 xmax: 1067494 ymax: 272710.9
CRS:            2263
ggplot() +
  geom_sf(data=zipcode_sf, aes(fill = borough)) 

# mapview function
# mapview(zipcode_sf)
mapview(altzip_sf)
# try leaflet with this
pal <- colorFactor(rainbow(5), zipcode_sf$borough)

leaflet(zipcode_sf) %>%
  addProviderTiles(providers$CartoDB.Positron) %>%  #add simpler basemap 
  addPolygons(stroke = FALSE, 
              smoothFactor = 0.3, 
              fillOpacity = 0.2,
    highlight = highlightOptions(color = "white", weight = 10, bringToFront = TRUE))  %>%
  addPolylines(color = "#444444", 
               weight = 1, 
               smoothFactor = 0.5, 
               opacity = 1.0, 
               highlightOptions = highlightOptions(color = "white", 
                                                   weight = 5, 
                                                   bringToFront = F, 
                                                   opacity = 1))
### combine zipcode geometric data with cleaned dataframe

#cleaned_sf <- geojsonsf::geojson_sf("combined_covid_data.geojson")

#colnames(zipcode_sf)[which(names(zipcode_sf) == "postalcode")] <- "ZIPCODE"

# make certain of zipcode overlap between the two dataframes
setdiff(cleaned_sf$ZIPCODE, altzip_sf$ZIPCODE)
character(0)
# altzip only need ZIPCODE and geometry columns:
altzip_sf_compact <- subset(altzip_sf, select = c('ZIPCODE', 'geometry'))

combined_sf <- merge(cleaned_sf, altzip_sf_compact, by = 'ZIPCODE')
Error in merge.sf(cleaned_sf, altzip_sf_compact, by = "ZIPCODE") : 
  merge on two sf objects not supported
# 
cleaned_sf_compact <- within(cleaned_sf, rm('geometry'))
cleaned_sf_compact <- cleaned_sf_compact %>% st_set_geometry(NULL)
class(cleaned_sf_compact)
[1] "data.frame"
#combined_sf <- merge(cleaned_sf, altzip_sf_compact, by = 'ZIPCODE')
#combined_sf <- merge(altzip_sf_compact, cleaned_sf_compact, by = 'ZIPCODE')

class(combined_sf)
[1] "sf"         "data.frame"
mapview(combined_sf)
# let's try leaflet with this new combined_sf:

pal <- colorFactor(rainbow(5), combined_sf$PO_NAME)

leaflet(combined_sf) %>%
  addProviderTiles(providers$CartoDB.Positron) %>%  #add simpler basemap 
  addPolygons(stroke = FALSE, 
              smoothFactor = 0.3, 
              fillOpacity = 0.2,
    highlight = highlightOptions(color = "white", weight = 10, bringToFront = TRUE))
sf layer is not long-lat datasf layer has inconsistent datum (+proj=lcc +lat_1=41.03333333333333 +lat_2=40.66666666666666 +lat_0=40.16666666666666 +lon_0=-74 +x_0=300000.0000000001 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs ).
Need '+proj=longlat +datum=WGS84'

#  addPolylines(color = "#444444", 
#               weight = 1, 
#               smoothFactor = 0.5, 
#               opacity = 1.0, 
#               highlightOptions = highlightOptions(color = "white", 
#                                                   weight = 5, 
#                                                   bringToFront = F, 
#                                                   opacity = 1))
LS0tCnRpdGxlOiBjb3ZpZDE5IFIgU2hpbnkgbm90ZWJvb2sKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKYGBge3J9CiMgbG9hZCBkZXBlbmRlbmNpZXMKbGlicmFyeShyZ2RhbCkKbGlicmFyeShnZW9qc29uaW8pCmxpYnJhcnkoZ2dwbG90MikKbGlicmFyeShsZWFmbGV0KQpsaWJyYXJ5KHBsb3RseSkKbGlicmFyeShSQ29sb3JCcmV3ZXIpCmxpYnJhcnkoZHBseXIpCmxpYnJhcnkoYnJvb20pCmxpYnJhcnkobWFwcHJvaikKbGlicmFyeShyZXNoYXBlMikKbGlicmFyeShkaWdlc3QpCmxpYnJhcnkobHVicmlkYXRlKQpsaWJyYXJ5KHNwZHBseXIpCmxpYnJhcnkocm1hcHNoYXBlcikKbGlicmFyeShyYXN0ZXIpCmxpYnJhcnkobWFwdmlldykKbGlicmFyeShzZikKYGBgCgoKYGBge3J9CiMgdHJ5IHRoaXMgd2l0aDogaHR0cHM6Ly93d3cuci1ncmFwaC1nYWxsZXJ5LmNvbS8xNjgtbG9hZC1hLXNoYXBlLWZpbGUtaW50by1yLmh0bWwKCnppcGNvZGVfYm91bmRhcmllcyA8LSByZWFkT0dSKCIuLi9aSVBfQ09ERV8wNDAxMTQvWklQX0NPREVfMDQwMTE0LnNocCIsCiAgbGF5ZXIgPSAiWklQX0NPREVfMDQwMTE0IiwgdmVyYm9zZSA9IEZBTFNFKQoKIyB0cnkgZGlmZmVyZW50IHRhYmxlIG9mIHppcGNvZGVzOgojYWx0X3ppcGNvZGVfYm91bmRhcmllcyA8LSByZWFkT0dSKCIuLi9hbHRfbnljemlwL255Y196aXBfY29kZV90YWJ1bGF0aW9uX2FyZWFzX3BvbHlnb25zLnNocCIsCiMgIGxheWVyID0gIm55Y196aXBfY29kZV90YWJ1bGF0aW9uX2FyZWFzX3BvbHlnb25zIiwgdmVyYm9zZSA9IEZBTFNFKQoKIyBmb3J0aWZ5IHRoZSBkYXRhIHRvIG9idGFpbiBkYXRhZnJhbWUgZm9yIHVzZSBieSBnZ3Bsb3QyCnppcGNvZGVfYm91bmRhcmllc19mb3J0aWZpZWQgPC0gdGlkeSh6aXBjb2RlX2JvdW5kYXJpZXMsIHJlZ2lvbiA9ICJaSVBDT0RFIikKCiMgcGxvdCBpdApnZ3Bsb3QoKSArCiAgZ2VvbV9wb2x5Z29uKGRhdGEgPSB6aXBjb2RlX2JvdW5kYXJpZXNfZm9ydGlmaWVkLAogICAgICAgICAgICAgICBhZXMoeCA9IGxvbmcsIHkgPSBsYXQsIGdyb3VwID0gZ3JvdXApLAogICAgICAgICAgICAgICBmaWxsID0gIiM2OWIzYTIiLCBjb2xvcj0id2hpdGUiKSArCiAgdGhlbWVfdm9pZCgpCmBgYAoKYGBge3J9CiMgdHJ5IHVzaW5nIGluZm8gZnJvbTogaHR0cHM6Ly9uY2Vhcy5naXRodWIuaW8vb3NzLWxlc3NvbnMvc3BhdGlhbC1kYXRhLWdpcy1sYXcvMy1tb24taW50cm8tZ2lzLWluLXIuaHRtbAoKY2xhc3MoemlwY29kZV9ib3VuZGFyaWVzKQoKZXh0ZW50KHppcGNvZGVfYm91bmRhcmllcykKCnBsb3QoemlwY29kZV9ib3VuZGFyaWVzLCBtYWluID0gIk5ZQyB6aXBjb2RlIGJvdW5kYXJpZXMiKQpgYGAKCgpgYGB7cn0KIyBjb252ZXJ0IHNwYXRpYWwgZGF0YSBmcmFtZSB0byBnZW9qc29uCiMgbnljX2pzb24gPC0gZ2VvanNvbl9qc29uKHppcGNvZGVfYm91bmRhcmllcykKIyBzYXZlIHRvIGxvY2FsIGZpbGUKIyBnZW9qc29uX3dyaXRlKG55Y19qc29uLCBmaWxlID0gIi4vbnljLmdlb2pzb24iKQoKYWx0X255Y19qc29uIDwtIGdlb2pzb25fanNvbihhbHRfemlwY29kZV9ib3VuZGFyaWVzKQpnZW9qc29uX3dyaXRlKGFsdF9ueWNfanNvbiwgZmlsZSA9ICIuL2FsdF9ueWMuZ2VvanNvbiIpCmBgYAoKYGBge3J9CiMgY3JlYXRlIGxlYWZsZXQgb2JqZWN0CnBhbCA8LSBjb2xvck51bWVyaWMoInZpcmlkaXMiLCBOVUxMKQoKbGVhZmxldChhbHRfemlwY29kZV9ib3VuZGFyaWVzKSAlPiUKIyAgYWRkUHJvdmlkZXJUaWxlcyhwcm92aWRlcnMkT3BlblN0cmVldE1hcCkgJT4lIAogIGFkZFRpbGVzKCkgJT4lICNhZGQgYmFzZW1hcAojICBhZGRQb2x5Z29ucyhzdHJva2UgPSBGQUxTRSwgc21vb3RoRmFjdG9yID0gMC4zLCBmaWxsT3BhY2l0eSA9IDEsCiMgICAgZmlsbENvbG9yID0gfnBhbChib3JvdWdoKSwKIyAgICBsYWJlbCA9IGJvcm91Z2gpICU+JQogIGFkZFBvbHlsaW5lcyhjb2xvciA9ICIjNDQ0NDQ0IiwgCiAgICAgICAgICAgICAgIHdlaWdodCA9IDEsIAogICAgICAgICAgICAgICBzbW9vdGhGYWN0b3IgPSAwLjUsIAogICAgICAgICAgICAgICBvcGFjaXR5ID0gMS4wKQpgYGAKCgpgYGB7cn0KIyMjIHJ1biBpdCBhbGwgYWdhaW4gYnV0IHRoaXMgdGltZSB3aXRoIFNGIChmYXN0ZXI/KQoKemlwY29kZV9zZiA8LSBzdF9yZWFkKCJhbHRfbnljemlwL255Y196aXBfY29kZV90YWJ1bGF0aW9uX2FyZWFzX3BvbHlnb25zLnNocCIpCmFsdHppcF9zZiA8LSBzdF9yZWFkKCIuLi9aSVBfQ09ERV8wNDAxMTQvWklQX0NPREVfMDQwMTE0LnNocCIpCgpnZ3Bsb3QoKSArCiAgZ2VvbV9zZihkYXRhPXppcGNvZGVfc2YsIGFlcyhmaWxsID0gYm9yb3VnaCkpIApgYGAKCmBgYHtyfQojIG1hcHZpZXcgZnVuY3Rpb24KIyBtYXB2aWV3KHppcGNvZGVfc2YpCm1hcHZpZXcoYWx0emlwX3NmKQpgYGAKCmBgYHtyfQojIHRyeSBsZWFmbGV0IHdpdGggdGhpcwpwYWwgPC0gY29sb3JGYWN0b3IocmFpbmJvdyg1KSwgemlwY29kZV9zZiRib3JvdWdoKQoKbGVhZmxldCh6aXBjb2RlX3NmKSAlPiUKICBhZGRQcm92aWRlclRpbGVzKHByb3ZpZGVycyRDYXJ0b0RCLlBvc2l0cm9uKSAlPiUgICNhZGQgc2ltcGxlciBiYXNlbWFwIAogIGFkZFBvbHlnb25zKHN0cm9rZSA9IEZBTFNFLCAKICAgICAgICAgICAgICBzbW9vdGhGYWN0b3IgPSAwLjMsIAogICAgICAgICAgICAgIGZpbGxPcGFjaXR5ID0gMC4yLAogICAgaGlnaGxpZ2h0ID0gaGlnaGxpZ2h0T3B0aW9ucyhjb2xvciA9ICJ3aGl0ZSIsIHdlaWdodCA9IDEwLCBicmluZ1RvRnJvbnQgPSBUUlVFKSkgICU+JQogIGFkZFBvbHlsaW5lcyhjb2xvciA9ICIjNDQ0NDQ0IiwgCiAgICAgICAgICAgICAgIHdlaWdodCA9IDEsIAogICAgICAgICAgICAgICBzbW9vdGhGYWN0b3IgPSAwLjUsIAogICAgICAgICAgICAgICBvcGFjaXR5ID0gMS4wLCAKICAgICAgICAgICAgICAgaGlnaGxpZ2h0T3B0aW9ucyA9IGhpZ2hsaWdodE9wdGlvbnMoY29sb3IgPSAid2hpdGUiLCAKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgd2VpZ2h0ID0gNSwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIGJyaW5nVG9Gcm9udCA9IEYsIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBvcGFjaXR5ID0gMSkpCmBgYAoKYGBge3J9CiMjIyBjb21iaW5lIHppcGNvZGUgZ2VvbWV0cmljIGRhdGEgd2l0aCBjbGVhbmVkIGRhdGFmcmFtZQoKI2NsZWFuZWRfc2YgPC0gZ2VvanNvbnNmOjpnZW9qc29uX3NmKCJjb21iaW5lZF9jb3ZpZF9kYXRhLmdlb2pzb24iKQoKI2NvbG5hbWVzKHppcGNvZGVfc2YpW3doaWNoKG5hbWVzKHppcGNvZGVfc2YpID09ICJwb3N0YWxjb2RlIildIDwtICJaSVBDT0RFIgoKIyBtYWtlIGNlcnRhaW4gb2YgemlwY29kZSBvdmVybGFwIGJldHdlZW4gdGhlIHR3byBkYXRhZnJhbWVzCnNldGRpZmYoY2xlYW5lZF9zZiRaSVBDT0RFLCBhbHR6aXBfc2YkWklQQ09ERSkKCiMgYWx0emlwIG9ubHkgbmVlZCBaSVBDT0RFIGFuZCBnZW9tZXRyeSBjb2x1bW5zOgphbHR6aXBfc2ZfY29tcGFjdCA8LSBzdWJzZXQoYWx0emlwX3NmLCBzZWxlY3QgPSBjKCdaSVBDT0RFJywgJ2dlb21ldHJ5JykpCgpgYGAKCgpgYGB7cn0KIyAKY2xlYW5lZF9zZl9jb21wYWN0IDwtIHdpdGhpbihjbGVhbmVkX3NmLCBybSgnZ2VvbWV0cnknKSkKY2xlYW5lZF9zZl9jb21wYWN0IDwtIGNsZWFuZWRfc2ZfY29tcGFjdCAlPiUgc3Rfc2V0X2dlb21ldHJ5KE5VTEwpCmNsYXNzKGNsZWFuZWRfc2ZfY29tcGFjdCkKI2NvbWJpbmVkX3NmIDwtIG1lcmdlKGNsZWFuZWRfc2YsIGFsdHppcF9zZl9jb21wYWN0LCBieSA9ICdaSVBDT0RFJykKYGBgCgpgYGB7cn0KI2NvbWJpbmVkX3NmIDwtIG1lcmdlKGFsdHppcF9zZl9jb21wYWN0LCBjbGVhbmVkX3NmX2NvbXBhY3QsIGJ5ID0gJ1pJUENPREUnKQoKY2xhc3MoY29tYmluZWRfc2YpCgptYXB2aWV3KGNvbWJpbmVkX3NmKQpgYGAKCmBgYHtyfQojIGxldCdzIHRyeSBsZWFmbGV0IHdpdGggdGhpcyBuZXcgY29tYmluZWRfc2Y6CgpwYWwgPC0gY29sb3JGYWN0b3IocmFpbmJvdyg1KSwgY29tYmluZWRfc2YkUE9fTkFNRSkKCmxlYWZsZXQoY29tYmluZWRfc2YpICU+JQogIGFkZFByb3ZpZGVyVGlsZXMocHJvdmlkZXJzJENhcnRvREIuUG9zaXRyb24pICU+JSAgI2FkZCBzaW1wbGVyIGJhc2VtYXAgCiAgYWRkUG9seWdvbnMoc3Ryb2tlID0gRkFMU0UsIAogICAgICAgICAgICAgIHNtb290aEZhY3RvciA9IDAuMywgCiAgICAgICAgICAgICAgZmlsbE9wYWNpdHkgPSAwLjIsCiAgICBoaWdobGlnaHQgPSBoaWdobGlnaHRPcHRpb25zKGNvbG9yID0gIndoaXRlIiwgd2VpZ2h0ID0gMTAsIGJyaW5nVG9Gcm9udCA9IFRSVUUpKQojICBhZGRQb2x5bGluZXMoY29sb3IgPSAiIzQ0NDQ0NCIsIAojICAgICAgICAgICAgICAgd2VpZ2h0ID0gMSwgCiMgICAgICAgICAgICAgICBzbW9vdGhGYWN0b3IgPSAwLjUsIAojICAgICAgICAgICAgICAgb3BhY2l0eSA9IDEuMCwgCiMgICAgICAgICAgICAgICBoaWdobGlnaHRPcHRpb25zID0gaGlnaGxpZ2h0T3B0aW9ucyhjb2xvciA9ICJ3aGl0ZSIsIAojICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgd2VpZ2h0ID0gNSwgCiMgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBicmluZ1RvRnJvbnQgPSBGLCAKIyAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIG9wYWNpdHkgPSAxKSkKYGBgCgo=